###################################
###Conjoint Analysis: Main#########
###################################


###################################
#Clean Environment
rm(list=ls())

#Load Libraries (Install if required)
pacman::p_load(cjoint, tidyverse, extrafont, haven, cregg, cairo)

#Import Fonts
font_import()
loadfonts(device="win")

df <- read_dta("conjoint_1_pooled.dta")



#Set up Data for Analysis
#####################################
#Set value labels (Marginal Means needs them to be specified)
df <- df %>%
  mutate(
    Early_retirement = factor(att1_earlyret,
           levels = c(1, 2, 3),
           labels = c("Expand early retirement for all", 
   "Leave early retirement unchanged", 
   "Abolish support for early retirement")),
    Childcare = factor(att2_childcare,
    levels = c(1, 2),
    labels = c("Expand childcare", "Leave childcare unchanged")),
    Inheritance_tax = factor(att3_inheritax,
          levels = c(1, 2, 3),
          labels = c("Increase inheritance taxes", 
  "Leave inheritance taxes unchanged", "Reduce inheritance taxes")),
    Immigration = factor(att4_immig,
      levels = c(1, 2, 3),
      labels = c("Regulated immigration with no upper limit", 
                 "Regulated immigration with upper limit", 
                 "Strong reduction of immigration")),
    Headscarf = factor(att5_headscarf,
    levels = c(1, 2),
    labels = c("Headscarf ban", "No headscarf ban")),
    Gender_quota = factor(att6_gender,
       levels = c(1, 2, 3),
       labels = c("50% gender quota", "30% gender quota", "No gender quota")),
    Ecology = factor(att7_ecology,
  levels = c(1, 2, 3),
  labels = c("Increase CO2-taxes strongly", 
             "Increase CO2-taxes moderately", "No increase of CO2-taxes")),
    Job_protection = factor(att8_jobprotection,
         levels = c(1, 2),
         labels = c("Expand job protection", "Leave job protection unchanged")),
    Housing = factor(att9_housing,
  levels = c(1, 2, 3),
  labels = c("Ban rent increases", "Slow down rent increases",
             "Leave rents unchanged"))
  ) %>%
  mutate(
    Early_retirement = factor(Early_retirement, 
           levels = c("Abolish support for early retirement", "Leave early retirement unchanged", "Expand early retirement for all")),
    Childcare = factor(Childcare, 
    levels = c("Leave childcare unchanged", "Expand childcare")),
    Inheritance_tax = factor(Inheritance_tax, 
          levels = c("Reduce inheritance taxes", "Leave inheritance taxes unchanged", "Increase inheritance taxes")),
    Immigration = factor(Immigration, 
      levels = c("Strong reduction of immigration", "Regulated immigration with upper limit", "Regulated immigration with no upper limit")),
    Headscarf = factor(Headscarf, 
    levels = c("Headscarf ban", "No headscarf ban")),
    Gender_quota = factor(Gender_quota, 
       levels = c("No gender quota", "30% gender quota", "50% gender quota")),
    Ecology = factor(Ecology, 
  levels = c("No increase of CO2-taxes", "Increase CO2-taxes moderately", "Increase CO2-taxes strongly")),
    Job_protection = factor(Job_protection, 
         levels = c("Leave job protection unchanged", "Expand job protection")),
    Housing = factor(Housing, 
  levels = c("Leave rents unchanged", "Slow down rent increases", "Ban rent increases"))
  )


##Interaction Variables: Relabel and factorize
# Relabel Oesch Classes
df <- df %>%
  mutate(oesch = case_when(
    class8_r == 1 ~ "Large employers & self-employed professionals",
    class8_r == 2 ~ "Small business owners",
    class8_r == 3 ~ "Technical professionals",
    class8_r == 4 ~ "Production workers",
    class8_r == 5 ~ "(Associate) managers",
    class8_r == 6 ~ "Office Clerks",
    class8_r == 7 ~ "Socio-Cultural Professionals",
    class8_r == 8 ~ "Service Workers",
    TRUE ~ NA_character_
  )) %>%
  mutate(oesch = factor(oesch, 
     levels = c("Large employers & self-employed professionals", 
                "Small business owners", "Technical professionals", 
                "Production workers", "(Associate) managers", 
                "Office Clerks", "Socio-Cultural Professionals", 
                "Service Workers")))

# Create Oesch Variable with EMC and PW
df <- df %>%
  mutate(oesch_2 = case_when(
    class8_r == 3 & educ_cat == 3 ~ "Educated Middle Class",
    class8_r == 7 & educ_cat == 3 ~ "Educated Middle Class",
    class8_r == 4 ~ "Production Workers",
    TRUE ~ NA_character_
  )) %>%
  mutate(oesch_2 = as.factor(oesch_2))

# Create Oesch Variable with SCP and PW
df <- df %>%
  mutate(oesch_3 = case_when(
    class8_r == 7 ~ "SCP",
    class8_r == 4 ~ "Production Workers",
    TRUE ~ NA_character_
  )) %>%
  mutate(oesch_3 = factor(oesch_3, levels = c("SCP", "Production Workers")))

# Create Oesch Variable with ESCP and PW
df <- df %>%
  mutate(oesch_4 = case_when(
    class8_r == 7 & educ_cat == 3 ~ "Educated SCP",
    class8_r == 4 ~ "Production Workers",
    TRUE ~ NA_character_
  )) %>%
  mutate(oesch_4 = as.factor(oesch_4))

# Create Age Group Variable with only two Groups
df <- df %>%
  mutate(age_groups = case_when(
    age <= 30 ~ "younger than 30",
    age > 59 ~ "60+",
    TRUE ~ NA_character_
  )) %>%
  mutate(age_groups = factor(age_groups, levels = c("younger than 30", "60+"))) %>%
  
  # Create Age Group Variable with four groups
  mutate(age_groups = case_when(
    age < 35 ~ "<35",
    age >= 35 & age <= 50 ~ "35-50",
    age > 50 & age <= 65 ~ "50-65",
    age > 65 ~ "65+",
    TRUE ~ NA_character_
  )) %>%
  mutate(age_groups = as.factor(age_groups))

# Create Sex Variable
df <- df %>%
  mutate(sex = case_when(
    female == 0 ~ "male",
    female == 1 ~ "female",
    TRUE ~ NA_character_
  )) %>%
  mutate(sex = factor(sex, levels = c("male", "female")))

# Create Education Variable
df <- df %>%
  mutate(educ_2 = case_when(
    educ_cat == 3 ~ "High",
    educ_cat == 2 ~ "Middle",
    educ_cat == 1 ~ "Low",
    TRUE ~ NA_character_
  )) %>%
  mutate(educ_2 = factor(educ_2, levels = c("Low", "Middle", "High")))

# Create Union Variable
df <- df %>%
  mutate(union = case_when(
    tradeunion == 3 ~ "No",
    tradeunion == 2 ~ "No",
    tradeunion == 1 ~ "Current Member",
    TRUE ~ NA_character_
  )) %>%
  mutate(union = factor(union, levels = c("No", "Current Member")))

# Create Urban Variable
df <- df %>%
  mutate(urban = case_when(
    urbanrural == 5 ~ "Countryside",
    urbanrural == 4 ~ "Village",
    urbanrural == 3 ~ "Small City",
    urbanrural == 2 ~ "Suburbs",
    urbanrural == 1 ~ "Big City",
    TRUE ~ NA_character_
  )) %>%
  mutate(urban = factor(urban, levels = c("Big City", "Suburbs", 
    "Small City", "Village", "Countryside")))

# Create Income Variable
df <- df %>%
  mutate(inc = case_when(
    income %in% 1:3 ~ "Low",
    income %in% 4:7 ~ "Middle",
    income %in% 8:10 ~ "High",
    TRUE ~ NA_character_
  )) %>%
  mutate(inc = factor(inc, levels = c("Low", "Middle", "High")))

# Make sure all previously coded variables are factors
df <- df %>%
  mutate(across(c(WC_male, WC, MC, young30, young35, potSD, vote_sd, young40,
                  Mage, country), as.factor))


# Filter out potential social democratic voters robust 1
df3 <- df %>% filter(ptv_sd > 5)

# Filter out potential social democratic voters robust 2
df <- df %>% mutate(ptvSD = ifelse(ptv_sd > 6 | lrs < 5, 1, 0))
df4 <- df %>% filter(ptvSD == 1)

# Filter out potential social democratic voters robust 3
df5 <- df %>% filter(lrs < 5)

# Four Groups according to Typology
df <- df %>% mutate(potentialgroups = case_when(
  ptv_sd > 5 & lrs < 5 ~ "Core supporters",
  ptv_sd > 5 & lrs > 4 ~ "Non-left Supporters",
  ptv_sd < 6 & lrs < 5 ~ "Left Non-Supporters",
  ptv_sd < 6 & lrs > 4 ~ "Non-Supporters",
  TRUE ~ NA_character_
))

#####################################################################################################
# Filter out potential social democratic voters
df1 <- df %>% filter(potSD == 1)

# Create Country Datasets (only potentials)
AT <- df1 %>% filter(country == "AT")
CH <- df1 %>% filter(country == "CH")
DE <- df1 %>% filter(country == "DE")
DK <- df1 %>% filter(country == "DK")
ES <- df1 %>% filter(country == "ES")
SE <- df1 %>% filter(country == "SE")

# Create Country Exclusion Datasets (only potentials)
No_AT <- df1 %>% filter(country %in% c("CH", "DE", "DK", "ES", "SE"))
No_CH <- df1 %>% filter(country %in% c("AT", "DE", "DK", "ES", "SE"))
No_DE <- df1 %>% filter(country %in% c("AT", "CH", "DK", "ES", "SE"))
No_DK <- df1 %>% filter(country %in% c("AT", "CH", "DE", "ES", "SE"))
No_ES <- df1 %>% filter(country %in% c("AT", "CH", "DE", "DK", "SE"))
No_SE <- df1 %>% filter(country %in% c("AT", "CH", "DE", "DK", "ES"))



####Conjoint Plot Adjustment
# Add H-Line to plot
hlines_values <- c(0, 4, 7, 11, 15, 18, 22, 26, 29, 33)



############################################################################
#Pooled Marginal Means only potentials (No interaction)
mm1 <- cj(df1, choice ~ Immigration + Headscarf + Gender_quota + Ecology + Childcare +
            Early_retirement + Inheritance_tax + Job_protection + Housing + country,  id = ~respondentid, estimate = "mm")

#Plot
cj1 <- plot(mm1, vline = 0.5, size = 5) +
  geom_point(position = position_dodge(0.75), color = "black", fill = "black") +  
  scale_x_continuous(limits = c(0.4, 0.6), breaks = c(0.3, 0.4, 0.5, 0.6, 0.7)) +
  theme_bw() + 
  xlab("Marginal Means") +    scale_color_grey(start = 0.1, end = 0.1, na.value = "red", aesthetics = "colour"   ) + 
  scale_color_grey(start = 0.1,
                   end = 0.1,
                   na.value = "red",
                   aesthetics = "colour"
  ) +
  theme_minimal() + 
  theme(text = element_text(family = "Times New Roman", size = 22)) + 
  theme(text = element_text(face = c('plain', 'plain','plain',
                  'bold','plain', 'plain', 'bold', 
                  'plain', 'plain', 'plain','bold'))) + 
  guides(fill = guide_legend(override.aes = list(color = NA)),
         color = FALSE,
         shape = FALSE) + 
  scale_y_discrete(labels = function(x) gsub("_", " ", gsub("\\(|\\)", "", x)))

#Remove NAs and Country Effects
cj1$data <- cj1$data[c(1:33),]
cj1$data <- na.omit(cj1$data) 

cj1 <- cj1 + geom_hline(yintercept = hlines_values, 
     color = "black", linetype = "dashed")
cj1

ggsave("Figure_5.tif", plot = cj1, device = "tiff", width = 16, height = 10)                    


############################
#marginal means potential vs. non potentials
mm_pot <- cj(df, choice ~ Immigration + Headscarf + Gender_quota + Ecology + Childcare +
            Early_retirement + Inheritance_tax + Job_protection + Housing,  id = ~respondentid, estimate = "mm", by = ~potSD)

#Plot
cj_pot_SD<-plot(mm_pot,group = "potSD", vline = 0.5, size=5) +
  geom_point(position = position_dodge(0.75)) +
  scale_x_continuous(limits = c(0.4, 0.6), breaks = c(0.3, 0.4, 0.5, 0.6, 0.7)) +
  theme_bw() + xlab("Marginal Means") +    scale_color_grey(start = 0.1, end = 0.1, na.value = "red", aesthetics = "colour"   ) + theme_minimal(base_size =
                  15) + theme(text=element_text(family="Times New Roman", size=22))+
  guides(color=guide_legend("Potential Social Democratic Voters")) + theme(legend.position="bottom") +
  theme(text=element_text(face = c('plain', 'plain','plain','bold',
                'plain', 'plain','bold',
                'plain', 'plain','plain','bold',
                'plain', 'plain','plain','bold',
                'plain', 'plain', 'bold',
                'plain', 'plain', 'plain','bold')))+ 
  geom_hline(yintercept = hlines_values, color = "black", linetype = "dashed") +
  scale_y_discrete(labels = function(x) gsub("_", " ", gsub("\\(|\\)", "", x)))
cj_pot_SD <- cj_pot_SD + aes(shape = factor(potSD))
cj_pot_SD <- cj_pot_SD + guides(colour = FALSE, na.translate = F) + theme(legend.position = "bottom") +
  scale_shape_discrete("", na.translate = F) +
  theme(text=element_text(face = c('plain')))
cj_pot_SD$data<-na.omit(cj_pot_SD$data)

ggsave("Figure_6.tif", plot = cj_pot_SD, device = "tiff", width = 16, height = 10)                    



#mm for educ_groups
mm_educ <- cj(df1, choice ~ Immigration + Headscarf + Gender_quota + Ecology + Childcare +
            Early_retirement + Inheritance_tax + Job_protection + Housing,  id = ~respondentid, estimate = "mm", by = ~educ_2)

#Plot
cj_educ<-plot(mm_educ,group = "educ_2", vline = 0.5, size=3) +
  geom_point(position = position_dodge(0.75)) +
  scale_x_continuous(limits = c(0.4, 0.6), breaks = c(0.3, 0.4, 0.5, 0.6, 0.7)) +
  theme_bw() + xlab("Marginal Means") +    scale_color_grey(start = 0.1, end = 0.1, na.value = "red", aesthetics = "colour"   ) + theme_minimal(base_size =
                  15) + theme(text=element_text(family="Times New Roman", size=22))+
  guides(color=guide_legend("Educational Groups")) + theme(legend.position="bottom") + geom_hline(yintercept = hlines_values, color = "black", linetype = "dashed") + 
  scale_y_discrete(labels = function(x) gsub("_", " ", gsub("\\(|\\)", "", x)))  +
  theme(text=element_text(face = c('plain', 'plain','plain','bold',
                'plain', 'plain','bold',
                'plain', 'plain','plain','bold',
                'plain', 'plain','plain','bold',
                'plain', 'plain', 'bold',
                'plain', 'plain', 'plain','bold')))
cj_educ <- cj_educ + aes(shape = factor(educ_2))
cj_educ <- cj_educ + guides(colour = FALSE, na.translate = F) + theme(legend.position = "bottom") +
  scale_shape_discrete("", na.translate = F) + 
  theme(text=element_text(face = c('plain')))

ggsave("Figure_7.tif", plot = cj_educ, device = "tiff", width = 16, height = 10)                    


#marginal means classes
mm_class <- cj(df1, choice ~ Immigration + Headscarf + Gender_quota + Ecology + Childcare +
            Early_retirement + Inheritance_tax + Job_protection + Housing,  
            id = ~respondentid, estimate = "mm", by = ~oesch)


#Plot
cj_class<-plot(mm_class,group = "oesch", vline = 0.5, size=3) +
  geom_point(position = position_dodge(0.75)) +
  scale_x_continuous(limits = c(0.38, 0.62), breaks = c(0.3, 0.4, 0.5, 0.6, 0.7)) +
  theme_bw() + xlab("Marginal Means") +    scale_color_grey(start = 0.1, end = 0.1, na.value = "red", aesthetics = "colour"   ) + theme_minimal(base_size =
                  15) + theme(text=element_text(family="Times New Roman", size=22))+
  guides(color=guide_legend("Class")) + theme(legend.position="bottom") +  scale_y_discrete(labels = function(x) gsub("_", " ", gsub("\\(|\\)", "", x))) +
  theme(text=element_text(face = c('plain', 'plain','plain','bold', 'plain', 'plain','bold', 'plain', 'plain','plain','bold',
                'plain', 'plain','plain','bold',
                'plain', 'plain', 'bold',
                'plain', 'plain', 'plain','bold')))


cj_class_2a <- cj_class + geom_hline(yintercept = hlines_values, color = "black", linetype = "dashed") + scale_y_discrete(labels = function(x) gsub("_", " ", gsub("\\(|\\)", "", x)))  +
  theme(text=element_text(face = c('plain', 'plain','plain','bold','plain', 'plain','bold','plain', 'plain','plain','bold','plain', 
                'plain','plain','bold','plain', 'plain', 'bold','plain', 'plain', 'plain','bold')))


#Only select relevant classes to display
cj_class_2a$data <- cj_class_2a$data[c(100:114, 133:147, 199:213, 115:132, 148:165, 214:231),]
cj_class_2a <- cj_class_2a + aes(shape = factor(oesch))
cj_class_2a <- cj_class_2a + guides(colour = FALSE, na.translate = F) + theme(legend.position = "bottom") +
  scale_shape_discrete("", na.translate = F) + 
  theme(text=element_text(face = c('plain')))

ggsave("Figure_8.tif", plot = cj_class_2a, device = "tiff", width = 16, height = 10)                    


#mm for age-groups
mm_age <- cj(df1, choice ~ Immigration + Headscarf + Gender_quota + Ecology + Childcare +
            Early_retirement + Inheritance_tax + Job_protection + Housing,  id = ~respondentid, estimate = "mm", by = ~age_groups)

#plot
cj_age<-plot(mm_age,group = "age_groups", vline = 0.5, size=3) +
  geom_point(position = position_dodge(0.75)) +
  scale_x_continuous(limits = c(0.4, 0.6), breaks = c(0.3, 0.4, 0.5, 0.6, 0.7)) +
  theme_bw() + xlab("Marginal Means") +    scale_color_grey(start = 0.1, end = 0.1, na.value = "red", aesthetics = "colour"   ) + theme_minimal(base_size =
                  15) + theme(text=element_text(family="Times New Roman", size=22))+
  guides(color=guide_legend("Age Groups")) + theme(legend.position="bottom") + geom_hline(yintercept = hlines_values, color = "black", linetype = "dashed") + 
  scale_y_discrete(labels = function(x) gsub("_", " ", gsub("\\(|\\)", "", x)))  + 
  theme(text=element_text(face = c('plain', 'plain','plain','bold',
                                   'plain', 'plain','bold',
                                   'plain', 'plain','plain','bold',
                                   'plain', 'plain','plain','bold',
                                   'plain', 'plain', 'bold',
                                   'plain', 'plain', 'plain','bold')))
cj_age <- cj_age + aes(shape = factor(age_groups))
cj_age <- cj_age + guides(colour = FALSE, na.translate = F) + theme(legend.position = "bottom") +
  scale_shape_discrete("", na.translate = F) + 
  theme(text=element_text(face = c('plain')))

ggsave("Figure_9.tif", plot = cj_age, device = "tiff", width = 16, height = 10)                    
                 